arXiv:1509.02894vl [physics.flu-dyn] 9 Sep 2015 


A General, Mass-Preserving Navier-Stokes Projection Method 


David Salac 

Mechanical and Aerospace Engineering, University at Buffalo SUNY, 318 Jarvis Hall, Buffalo, NY 14260-4400, USA 


Abstract 

The conservation of mass is common issue with multiphase fluid simulations. In this work a novel projection method 
is presented which conserves mass both locally and globally. The fluid pressure is augmented with a time-varying 
component which accounts for any global mass change. The resulting system of equations is solved using an efficient 
Schur-complement method. Using the proposed method four numerical examples are performed: the evolution of a static 
bubble, the rise of a bubble, the breakup of a thin fluid thread, and the extension of a droplet in shear flow. The method 
is capable of conserving the mass even in situations with morphological changes such as droplet breakup. 
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1. Introduction 

A wide variety of techniques have been proposed to 
model multiphase fluid systems with incompressible, im¬ 
miscible fluids. These include explicit front-tracking tech¬ 
niques ULi , the volume-of-fluid method ii , the ^ase- 
fleld method [I, 0 , and the level set method [1,0. A 

constant challenge in each of these techniques is the conser¬ 
vation of mass during the course of the simulation. Each of 
these methods handle this challenge differently. For exam¬ 
ple, the volume-of-fluid methods has excellent mass conser¬ 
vation properties 0 at the expense of requiring complex 
heuristic interface reconstructions techniques to calculate 
geometric quantities such as curvature [ijJ, ll2|, • 

Unfortunately, there are many physical systems where 
high accuracy of geometric quantities are required. One 
example are models of liposome vesicles where inter facial 
forces depend on high order derivatives of the interface’s 
curvature [ij, [ll, ll6|. Unlike volume-of-fluid techniques, 
front tracking, phase-field, and level-set methods are able 
to provide higher geometric accuracy. This accuracy is 
obtained at the expense of natural volume conservation 
and thus special care must be taken to ensure that mass 
does not change over the course of a simulation. 

There have been numerous attempts to improve the 
mass conservation of such methods. For example, level 
set methods have been adjusted by Lagrangian particles 
which are used to correct the level set function [l7|, ll8| or 
level sets have been combined with volume of fluid meth¬ 
ods [l0 . Other simply shift the interface function to match 
the volume constraint (2^ . 

The major issue with these types of corrections is that 
the interface becomes decoupled from the underlying flow 
field. As an example consider an interface on a uniform 
Cartesian grid. Fig. [H The interface begins on the left side 


and the underlying fluid flow-field dictates that it should 
move to the right one grid spacing. After this time step 
it is determined that errors in the simulation resulted in a 
mass gain. A simple and often-used correction is to simply 
move the interface to obtain mass conservation. In this 
case the effect is the interface will not move the full amount 
that the underlying flow field dictates. The movement of 
the interface and the underlying flow field have become 
decoupled. Any externally applied correction for front¬ 
tracking, phase-field or level set methods will demonstrate 
similar behavior. 



Figure 1: Schematic of velocity-interface decoupling. The interface 
begins on the left (solid line). The underlying flow fleld dictates 
that the interface move to the right and should end as the dash-line 
after a single time step. Due to an externally applied correction the 
interface instead ends the time step at the dotted line. 


In this manuscript a different approach is taken. Instead 
of adjusting the interface to achieve mass conservation a 
novel Navier-Stokes projection method is developed which 
ensures mass conservation. Unlike the previous method 
the interface is simply advected due to the underlying 
flow-field, it is the flow-field itself which explicitly takes 
into account any possible mass-loss. As will become ap- 
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parent later in the manuscript it is useful to think of this 
method as modifying the pressure so that it can handle 
not only local incompressibility but also global mass con¬ 
servation. Note that while this manuscript will focus on 
a particular Navier-Stokes numerical implementation, the 
concept presented here extends to any type of multiphase 
fluid simulation that uses a projection method. 

The remainder of the manuscript is as follows. In Sec¬ 
tion the single-fluid formulation of multiphase fluid flow 
is briefly described. The novel mass-preserving Navier- 
Stokes projection method is described in Section [3l The 
numerical implementation is given in Section 01 which is 
followed by numerical experiments in Section [5l A short 
conclusion is presented in Section [6l 

2. Single-Fluid Navier-Stokes Equations 

In this sections a brief introduction to the single-fluid 
formulation of multiphase fluid flow is presented. In a mul¬ 
tiphase fluid system the interface between two immiscible 
fluids evolves over time. In this work a level-set descrip¬ 
tion of the interface is used. Let the evolving interface 
be given as the set of points where the level-set function 
is zero, r(t) = {x : (t){x^t) = 0}. The evolution of the 
level set function (j){x,t) over time will implicitly deter¬ 
mine the location of the interface. Following convention 
the interior fluid, is given by 0 < 0 while the outer 
fluid, is given by 0 > 0, Fig. The entire domain is 
given as Using the level-set description it is 

possible to obtain geometric information of the interface 
easily. For example, the outward facing normal is simply 
n = V0/II V(/)|| while the total curvature can be calculated 
8iS K = \/ ■ n. 



Figure 2: The level-set description of a multiphase fluid system. The 
outward facing normal is shown for clarity. 

In each domain the Navier-Stokes equations hold, 

= -Vp± + V • (1) 

V • = 0, (2) 

where p is the fluid density, p is the fluid viscosity and g is 
any body force term, such as gravity. The fluid equations 


are coupled by a jump in the stress at the interface, 

[-pn + p (Vn + V^n)] n = f, (3) 

where [ ] indicates the jump of a quantity (outside minus 
inside) across the interface and / are any forces, such as 
tension, which act on the interface. 

A difficulty in multiphase fluid simulations is the solu¬ 
tion of this set of coupled but discontinuous differential 
equations. Some techniques, such as the Immersed Inter¬ 
face Method [111, augment the discretization of the dif¬ 
ferential equations to take into account the jumps across 
the interface. Another technique, which is explained here, 
is to model the domain as a “single” fluid with spatially 
varying properties Blii. The interface condition, Eq. 
m, is accounted for by converting the singular force con¬ 
tribution at the interface into a body-force term localized 
around the interface. 

Define the smooth Heaviside, function as 




1 + f + i sin 



(j) < —e 

\4>\ < £ 

(j)> e 


( 4 ) 


where 5 is proportional to the grid spacing. From the 
definition of the Heaviside function define the smoothed 
Dirac-Delta function as 5s{(j)) = dH^{(j)) jd(j)^ 


SM = 



l + cos(f) 


l<^l > e 

\(t>\<£' 


( 5 ) 


The use of the Heaviside and Dirac functions allows for 
the calculation of the volume enclosed by a given level set 
as well as the surface area as a integrals over the domain. 


V{^) = 

Jn 

( 6 ) 

A{4>) = 

[ SM\W4>\\ dV. 

Jn 

( 7 ) 


Using the Heaviside function the density and viscosity 
are given by Pe(0) = p~ — p~) H^{(j)) and p£{(j)) = 

p~ + (p+ — p~) These functions also allow for the 

transformation of singular interface forces into localized 
body force terms. For example, let the only singular inter¬ 
facial force be from a uniform surface tension, f = a/^n, 
where a is the coefficient of surface tension. In conjunc¬ 
tion with the smoothed density and viscosity definitions 
this results in a single Navier-Stokes equation valid in the 
entire domain, 

Du 

= - Vp + V ■ {pe{4>) (V-U + V'^u)) 

- aSei4>)KV4> + 


2 


V • u =0, 


(8) 

( 9 ) 









where the body force term has been written as a smooth 
function of the level-set. This type of single-fluid formula¬ 
tion for multiphase flow has been used to model bubbles 
and droplets 23, 24| and vesicles [mu. 

The results below focus on bubbles and droplets under 
the influence of surface tension and gravity. Therefore, the 
interface and body force are restricted to this particular 
case. The dimensionless form of the single-fluid Navier- 
Stokes equations can be obtained by defining a character¬ 
istic length /o, and velocity i^o, from which a characteristic 
time can be obtained, to = Iq/uq- The density and viscos¬ 
ity are normalized by the values of the outer domain. 


p{^) = X + 0-X)H{ct>), (10) 

+ (</.), (11) 


where A = p~ /p'^ and p = p~ /are the density and 
viscosity ratio, respectively, and the 5 has been dropped 
for clarity. 

Define the the Reynolds number as 


P'^IqUo 
Re =-—, 


the Weber number as 


We = 


P'^lpul 

(J 


and the Fronde number as 


Fr = 


^0 

^0^0 ’ 


( 12 ) 


(13) 


(14) 


where go is the strength of gravity pointing in the g direc¬ 
tion. The Navier-Stokes equations now become 


+ (15) 

Note that in this particular formulation gravity is written 
as (p(0) — 1)^ so that uniform acceleration of the entire 
domain does not occur. 


3. A Mass-Preserving Projection Method 

Projection methods are a common technique to solve 
the Navier-Stokes equations [IgI. In such schemes a tenta¬ 
tive velocity field is first determined and then the pressure 
is calculated to enforce the local volume conservation. In 
addition to conservation of local volume, Eq. m, it would 
be advantageous to explicitly enforce global volume con¬ 
servation. Even if local volume conservation is performed 
exactly errors in the advection of the interface can intro¬ 
duce unwanted volume change over the course of a sim¬ 
ulation. In this novel projection method both local and 
global volume conservation are enforced. 


The first step is to determine a tentative velocity field 
at time in the absence of the pressure and forces. In 
this work the semi-Lagrangian formulation is used: 

7/* _ qi'Tl 1 

p m 

The departure velocity, for a grid point x is determined 
by tracing characteristics backward in time, 

Xd = X — Atu^ (x ), (17) 

Ud = Pu{xd,t"'), (18) 


where Pu is an interpolant of the velocity field at time 
. Additional information regarding the semi-Lagrangian 
method for Navier-Stokes equations is given in Ref |27j . 

The next step of a standard projection method would 
be to determine the pressure to satisfy local volume con¬ 
servation, Eq. (|9]). In this work global conservation is also 
considered. Global volume conservation can be enforced 
by looking at the rate at which the enclosed volume, E, 


will change over time [15 


L 


n ■ u^+'^dA = ^ 
dt 


(19) 


for an outward facing unit normal, n, to the interface F. 
Unfortunately, the pressure field alone can not satisfy both 
constraints. Therefore the pressure is split into a constant 
and spatially-varying component: 


p = p+(1 - H{(t)))po. (20) 


The constant pressure component, po, only has a contri¬ 
bution to the overall pressure field in the region given by 
0 < 0 and can be thought of enforcing global volume con¬ 
servation. The spatially-varying component, p, takes the 
place of the standard pressure and will be used to enforce 
local volume conservation. Expanding the pressure com¬ 
ponent results in 


-Vp = -v(p + (1 - H{r))po) = -Vp + S{r)po^r, 

( 21 ) 

which demonstrates that po only has a contribution near 
the interface. 

Using this definition of pressure the projection step is 
now 

1 6{(t)^)poV(t)^ 

At p{4>^) {p) 

_ 5{r)i^v4> , pm-1 - .221 

{p)We p(</.«)Fr ^ ’ 

Eor quantities localized around the interface the average 
density is used, (p) = (p+ + P~)l‘^' 

To determine p and po both conservation conditions are 
applied to the projection step, Eq. (j22j). Applying local 
conservation, Eq. m, results in 


}u* .5(0")p(0")-l \ 

\At {p)We p{4>^)Ft ^) ■ 


(23) 
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Global conservation can be enforced by combining Eqs. 
m and ([2^ and also requiring that the time-evolution of 
the volume be such that any previous errors are corrected, 
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(24) 


where is the initial enclosed volume. 

Let p represent the vector holding the the values of p 
over the entire discretized domain. The surface integrals in 
Eq. (|M|) can be written as a summation over a discretized 
domain, J-^fdA ^ where grid points are 

given by = d{(j){xi^j^k)), fij,k is the function 

value at the grid points, and dV is the volume of each cell 
surrounding a grid point. It is then possible to define the 
following linear operators: 



» ip, 

(25) 

^V-(<5(0")V0). 

-Pol, 

(26) 

[ —^—n ’Vp dA ^ 

Jr p{4>r 

T ~ 

^ s p 

(27) 


- Poa. 

(28) 


The two constraint equations, Eqs. ([23]) and 
be written in block matrix-vector form. 


can now 


L 1 


P 


d 

T 

s a 


Po. 


e 


(29) 


where d is the discretization of the right-hand-side of Eq. 
(|23]) while e is the evaluation of the integral on the right- 
hand-side of Eq. (HH). This particular block matrix form 
can be solved efficiently through the use of a Schur decom¬ 
position. 


d 

(30) 

The Schur complement is a rank-1 update on the matrix 




s-^ 

0 


0 

^-1 


I -a-H 
0 1 


L: 

S = L-hs^, (31) 

a 

which has an inverse given by the Sherman-Morrison for¬ 
mula. 


S-'^ =L-'^ + 




(32) 


So long as a — L~^l is non-zero the inverse is defined. 
This condition is not proven here, but a check is performed 
during all simulations and has never been violated. Note 
that as the method as described does not depend on any 
particular spatial discretization. The particular discretiza¬ 
tion used in the numerical experiments below will be dis¬ 
cussed in Sec. 


4. Numerical Implementation 

In this section the particular numerical implementation 
used in the numerical experiments of Sec. [5] is briefly dis¬ 
cussed. In particular the description and advection of the 
interface and the discretization of the projection method 
is presented. 

4 . 1 . Interface Description and Advection 

The interface separating the multiple fluid phases is de¬ 
scribed here using a gradient-augmented level set method 
[n, 113 . In this extension of the level set method the 
gradient field of the level set is explicitly tracked in addi¬ 
tion to the level set function. This allows for the accurate 
determination of interfacial quantities such as curvature 
m, [ 3 ^ . In particular the results shown below use a re¬ 
cent semi-implicit extension of the original level set Jet 
scheme [sij. 

Consider the implicit tracking of an interface by way of 
a level set function: r(t) = {x : (j){xA) = 0}. In addition 
to the level set also track the gradient field of the level set, 
-0 = V0. This allows for the determination of Hermite 
interpolants using only information local to a cell on a 
Cartesian grid [29|. Unfortunately, the original Jet scheme 
does not work well for stiff advection equations, such as 
those involving surface tension, and the recent Semi Jet 
method was developed in response 0. 

The Semi Jet method begins by performing a semi- 
implicit, semi-Lagrangian update on the level set field, 

^ - /3V2<()", (33) 

where (j)^ is the level set value at the departure location 
and /3 = 0.5 is a constant. 

To update the gradient field note that the smoothing 
operation can be captured using a smoothing operator, 

- /3VV", (34) 

where and are both known at grid points after 

solving Eq. (issii . At each grid point define a sub-grid of 
spacing e ^ h, where h is the grid spacing: 

x^ = X ^ qe ioT q ^ { — 1,1}^, (35) 

with p = 2 in two-dimensions and p = 3 in three- 
dimensions. This results in a sub-grid of either 4 or 8 
points, depending on the dimension of the simulation. On 
each of these sub-grid points an updated level-set value 
is obtained using an explicit semi-Lagrangian update with 
the smoothing source term given in Eq. (El), 

= (36) 

Once the level set values on the subgrid are calculated 
updated gradient values at grid points can be obtained 
using finite difference approximations. Eor example, in 
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two-dimensions the gradient in the x-direction at a grid is 
given by _ 0 (-id) / 4 e. 

Curvature values are calculated at grid points using 


K = 



(37) 


in two-dimensions with a similar expression for three- 
dimensions. The values used to calculate the curvature 
were calculated by averaging derivatives. For example, 
when computing 4^xx and (pxy for use in the curvature 
calculation the following averaging would be used: 


= \{r + D,4>) , (38) 

4^xx = \ {Dx'ip^ + Dxx(t>) j (39) 

^xy = ^ {Dy^p^ -h Dx'lp^ + DxyCp) , (40) 


5. Numerical Experiments 

The mass-preserving projection method described above 
is tested using four different multiphase fluid systems. 
These are: 1) a static bubble, 2) a rising bubble, 3) a fluid 
thread undergoing a Plateau-Rayleigh instability, and 4) 
a droplet in shear flow. In all cases the mass-preserving 
method, denoted as the conserving scheme, is compared 
to a simulation without the method described above. This 
non-conserving scheme is computed by using the same dis¬ 
cretization as the conserving scheme with the pressure cal¬ 
culated as the solution to Eq. (|23|) with po = 0. All of the 
results will be three-dimensional results using uniform grid 
spacing, although some figures will show two-dimensional 
slices instead of three-dimensional figures. 

5.1. A Static Bubble 


where Dx^ Dy^ Dxx and Dxy are the finite difference ap¬ 
proximations for the derivatives and and are the 
components of the gradient field tracked by the Jet. The 
curvature was then extended away from the interface in the 
normal direction. Complete information about the Semi- 
Jet method and curvature calculation can be found in Ref. 


4 . 2 . Discretization of the Projection Method 

The Navier-Stokes equations and projection method are 
discretized on a collocated Cartesian mesh using stan¬ 
dard finite difference approximations. The variable co¬ 
efficient Poisson equation, Eq. (ESI, is discretized using 
compact second-order finite difference approximations. In 
one-dimension this is written as 


^ Pih/2 fe+1 - Pi) - Pi-1/2 jPi - P»-l) 

* h‘^ 


where the inverse density at half-grid points is obtained 
using harmonic averaging. 


-1 _ Pi A- Pi+1 

- 2p,p,+, ■ 


(42) 


It was found that the harmonic averaging improved the 
stability for large density variations. Similar expressions 
hold in other dimensions. The remaining discretizations in 
Eqs. (|2^ - (|28|) were also performed using centered second- 
order finite difference approximations. 

Note that the method described here is an approximate 
pressure discretization. It is well known that this partic¬ 
ular discretization does not suffer from the pressure de¬ 
coupling effect which can occur from a discrete discretiza¬ 
tion on a collocated mesh [n On the other hand 

the fact that the divergence-free velocity field condition 
is not discretely enforced means that small mass errors 
will accumulate over time. An approximate discretization 
was purposefully chosen to demonstrate the ability of the 
proposed method to conserve mass even if the underlying 
method is not discretely mass conserving. 


The first numerical experiment is the evolution of a 
droplet in the absence of gravity, Eig. [3l The viscosity 
ratio is set to 7 ^ = 0.01 while the density ratio is A = 0.001. 
The Reynolds number is 10 and the Weber number is 1. 
Gravity is ignored, which is equivalent to setting Er = oc. 
The initial shape is an ellipsoidal interface with an axis of 
1 in the x— and 2 :—directions and 2 in the ^—direction. 
The computational domain is a cube given by [—4,4]^ and 
a 65^ grid, resulting in a grid spacing of h = 0.125. The 
time step is At = 0.01. 




(a) t =0 


(b) t=0.5 


(c) t=l 





(d) t =2 


(e) t =10 


(f) t =20 


Figure 3: Evolution of a bubble under surface tension with Re = 10, 
We = 1, and Fr = 00 . The density ratio is A = 0.001 while the 
viscosity ratio is 77 = 0.01. The final sphere radius is 1.261. 


Under the influence of surface tension the bubble evolves 
into a sphere with the same volume. It is therefore ex¬ 
pected that the final shape will be a sphere with a radius 
of 1.26. The simulated final radius is 1.261, very close to 
the theoretical value. 

The distribution of the spatially-varying pressure com¬ 
ponent, p, and the total pressure p along the x—axis for 
y = 0 and 2 : = 0 at a time of t = 20 are shown in Eig. 
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[H Results are similar in the other directions. It is clear 
that at steady-state the majority of the pressure is due 
to P 07 which has a value of po = 1.584. According to the 
Laplace-Young condition it is also expected that the dif¬ 
ference between the inner and outer pressure should be 
p~ — = 2/^/We. Using an analytic final radius of 1.26 

and the simulation parameters this results in an expected 
pressure jump of 1.587, which matches well with the sim¬ 
ulated final pressure. 
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(b) Total Pressure: p = p {1 — H{(f)))pQ 


Figure 4: The pressure along the x-axis with y = 0 and 2 ; = 0 
at a time of t = 20 for the static bubble shown in Fig. [S] The 
spatially varying pressure component, p, is nearly zero while the 
pressure constant is po = 1.584. This is very near to the predicted 
static pressure of 1.587. 


A comparison between the mass conserving and non¬ 
mass conserving simulations is given in Fig. [5l Both 
the volume and volume-error over the course of the sim¬ 
ulation are presented. The volume error is computed as 
(U(t) —U( 0 ))/U( 0 ) where U(t) is the instantaneous volume 
and V (0) is the initial volume. The use of an approximate 
projection method, combined with the non-conservation 
properties of the level set method, results in excessive mass 
loss in the non-conserving simulation. This mass loss re¬ 
sults in the complete disappearance of the bubble in the 
non-conserving simulation, while the conserving scheme 
has a maximum volume error of approximately 2 x 10 “^. 



(a) Volume 



(b) Volume Error 

Figure 5: Comparison between volume preservation in the conserv¬ 
ing and non-conserving schemes for a static bubble with the initial 
condition shown in Fig. [3] After t = 15 the bubble has completely 
disappeared in the non-conserving simulation. 


5.^. A Rising Bubble 

The next numerical experiment is that of a rising bub¬ 
ble. The computational domain is a cube given by [—5, 5]^ 
using a 101^ grid, giving a grid spacing of 0.1. The time 
step is chosen to be 0.01. The initial bubble is a sphere 
centered at the origin with a radius of 1 . 

To allow for the long-time evolution of the bubble to be 
modeled after every time step the level set Jet and com¬ 
plete fluid field are re-centered in the computational do¬ 
main using a cubic interpolation procedure. Using the val¬ 
ues for the spherical cap Bubble A in Table 1 of Hnat and 
Buckmaster [sj the density ratio is set to A = 0.001142 
while the viscosity ratio is 77 = 0.008474. Based on the di¬ 
mension of the experimental bubble, the final rise velocity, 
and the given surface tension the dimensionless parameters 
are Re = 9.694, We = 7.6375 and Fr = 0.7754. 

As the bubble rises it evolves into the expected spherical 
cap, see Fig. [ 6 ] for a three-dimensional figure and Fig. [7| 
for the cross-section in the x — y plane. The experimen¬ 
tal results indicate that the aspect ratio of the bubble at 
steady-state should be 2.7, which compares to a value of 
2.8 from the simulation. With the given parameter set it 
is expected that the dimensionless rise velocity should be 
equal to one. The rise velocity over the course of the sim- 
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(d) t=4 


(e) t =10 


(f) t=20 


Figure 6 : Evolution of a three-dimensional bubble with density ratio 
of A = 0.001142 and viscosity ratio of 77 = 0.008474 with Re = 9.694, 
We = 7.6375 and Fr = 0.7754. These properties correspond to Bub¬ 
ble A in Table 1 in Hnat and Buckmaster (a. Experimental results 
show the aspect ratio should be 2.7 while the simulation results in 
an aspect ratio of 2 . 8 . 


O 


(a) t =0 (b) t=l (c) t =2 


Figure 8 : The rise velocity of the center of mass for the three- 
dimensional bubble in Figs. [ 6 ] and [7| The experimental results of 
Hnat and Buckmaster indicate the final dimensionless rise veloc¬ 
ity should be 1. The simulation results in a rise velocity of 0.9922. 
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Figure 7: Evolution of the x — y plane of a three-dimensional bubble 
for the result shown in Fig. ( 6 ] 


(a) Volume 



ulation is given in Fig. [HI with a steady-state rise velocity (b) Volume Error 

of 0.9922. 


The evolution of the volume over the course of the sim¬ 
ulation is presented in Fig. |9l Note that due to discretiza¬ 
tion errors the initial volume is not equal to 47r/3 ^ 4.19 
but is calculated numerically as 4.252. As with the static 
bubble example both the volume and the volume error 
are given. Over the course of the simulation the volume 
for the non-conserving scheme varies from a minimum of 
4.13 to a maximum of 4.39, while the conserving scheme 
only varies slightly. It should also be obvious that if the 
simulation was allowed to continue volume errors in the 
non-conserving scheme would continue to grow, which is 
not true for the conserving scheme. 


Figure 9: Comparison between volume preservation in the conserv¬ 
ing and non-conserving schemes for a rising bubble with the initial 
condition shown in Figs. [6]and[7| The volume gain will continue in 
the non-conserving scheme while the volume will remain constant in 
the conserving scheme. 


5.3. Plateau-Rayleigh Instability 

As a third test case consider the evolution of a fluid 
thread. It is well known that if the shape has an initial 
perturbation this thread will undergo a Plateau-Rayleigh 
instability and split into separate droplets [35[. In this 
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case the computational domain is a rectangle spanning 
[-5,5] X [-10,10] X [—5,5] while the number of grid points 
is 101 X 201 X 101 and the time step is At = 0.05. Periodic 
boundary conditions are given in all directions. 

Initially a thread is aligned along the ^-direction with 
an initial radius given by 1 + 0.2 cos(0.27r^) in the x — z 
plane. It is then allowed to evolve under the conditions 
of Re = 10, We = 1, and Fr = oo with a density ra¬ 
tio of A = 10 and viscosity ratio of r] = 10. The three- 
dimensional evolution of the interface is given in Fig. [TOl 
As time progresses the initially thin regions of the thread 
become thinner. At a time of t = 28 the thread is about 
to pinch off and extremely thin regions are observed be¬ 
tween individual droplets, which are visible at t = 28.5. 
Over time the droplets become spherical to minimize the 
surface tension. 



* 

-10 -2 

(c) t=28 (d) t=28.5 



Figure 10: Evolution of a fluid thread with a density ratio of A = 10, 
viscosity ratio of 77 = 10, Re = 10, We = 1, and Fr = 00 . Periodicity 
is assumed in all directions. With the given initial shape the droplet 
undergoes a Plateau-Rayleigh instability and breaks into spherical 
individual droplets. 


A comparison between the conserving and non¬ 
conserving schemes for the Plateau-Rayleigh instability is 
presented in Fig. [TTl As with the static bubble case the use 
of the non-conserving scheme results in complete volume 
loss. The conserving scheme has a maximum mass loss 
during the breakup observed at approximately t = 28.5, 
but quickly recovers to an error of less than 10“^. 



(a) Volume 



(b) Volume Error 

Figure 11: Comparison between volume preservation in the con¬ 
serving and non-conserving schemes for a fluid thread undergoing a 
Plateau-Rayleigh instability with the initial condition shown in Fig. 
CS] Due to volume loss the droplets in the non-conserving scheme 
disappear. 


5.4- Droplet in Shear Flow 

As a final test consider the behavior of a droplet in the 
presence of shear flow. Under certain situations the droplet 
will extend and “daughter” droplets may form at the tip of 
the original droplet [36|. Here the computational domain 
spans [—16,16] x [—2, 2] x [—2, 2] with a 513 x 65 x 65 grid. 
The time step is At = 0.025. Periodic boundary conditions 
are taken in the x— and 2;— directions while wall-boundary 
conditions are taken in the ^—direction. The shear flow is 
imposed by specifying a velocity of {y, 0, 0) on the y = 2 
and y = —2 walls. 

The result for matched density and viscosity, A = 1 and 
77 = 1, with Re = 1, We = 1, and Fr = oc for an initially 
spherical droplet of radius 1 is given in Fig. [T2j This 
figure shows the x — y cross-section of the droplet up to 
the formation of the first “daughter” droplets from the 
tips. 

Tracking the volume and volume error of the course of 
the simulation. Fig. [131 demonstrates results similar to the 
previous numerical tests. Over time the errors in the non¬ 
conserving scheme accumulate and results in the complete 
disappearance of the droplet. The conserving scheme, on 
the other hand, is capable of ensuring that the mass errors 
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Figure 12: Evolution of a droplet in the x — y plane under shear flow. 
The density ratio is A = 1, the viscosity ratio is 77 = 1, while Re = 1, 
We = 1, and Fr = 00 . The normalized shear rate is 1. 




(b) Volume Error 

Figure 13: Comparison between volume preservation in the conserv¬ 
ing and non-conserving schemes for the droplet in shear flow with the 
initial condition shown in Fig. [121 Due to volume loss the droplet in 
the non-conserving scheme disappears. 


do not grow over time. 

6. Conclusion 

In this work a new mass conserving projection method 
for multiphase Navier-Stokes systems has been presented. 
Various numerical tests have demonstrated that the pro¬ 
posed scheme is capable of conserving the mass in a wide 
variety situations, including a static bubble, a rising bub¬ 
ble, a fluid thread undergoing a Plateau-Rayleigh instabil¬ 
ity and a droplet in shear flow. Unlike many other mass 
correction schemes the material interface advects with the 
underlying fluid field. The general idea is applicable to 
any system being investigated using a projection-based 
Navier-Stokes solver and is not limited to the particular 
discretization used here. 
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